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THE METHOD OF AVERAGES APPLIED TO THE 


KS DIFFERENTIAL EQUATIONS 
by 

Otis F. Graf, Jr.. Alan 'Iu<*Ilfr and Stophc»n Siarko 


1.0 INTRODUCTION 


Hnlutions to ordinary di f 1 (^rontlal equations are usually 
obtained through the application of either analytical or nu- 
merical solution methods. Each approach has its own set of 
advantages and disadvantages. A particular method is usually 
chosen according to how it may satisfy foreseen applications. 

Analytical mf'thods require extensive development of math- 
c'matical formulas. Once they have been obtained, the formulas 
can be used to study certain global properties el the solu- 
tion. Also, the analytical formulas can be used in computing 
machines to provide extremely rapid numerical calculations. 

There is relatively little apriori mathematical develop- 
ment of numerical methods. They allow t'xtromely accurate* iiu- 
mericul calculations of the solution and are usually not dif- 
ficult to program on the computer. However, the solution must 
be developed step-by-step, which can lake large amount.s of 
computer time. Also, very littlt? cjualitative or global inft^r- 
mation on the solution is available. 

For some applications, it may be advantageous to devf’lop 
now specialized methods that combine features of both the ana- 
lytica. and the numerical methods. This is the suggested ap- 
proach for those applications that require repetitive solutions 
of the differential equations, but for which no analytical so- 
lutions are available. Therefore, the advantage in speed of 
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an analytical solution might be pax'tially realized by using 
a limited application of an analytical method. Numerical 
step-by-step calculation? still need to be done, but usually 
in a more efficient manner. 

This report concerns a new approach for the solution of 
artifical satellite trajectory problems. The basic idea is 
to apply an analytical solution method (the Method of Averages) 
to an appropriate formulation of the orbital mechanics equa- 
tions of motion (the KS-element differential equations). The 
result is a set of transformed equations of motion that are 
more amenable to numerical solution. 

The following three subsections in this Introduction 
give an overview of the Method of Averages, a review of the 
KS-elements, a short discussion of orbit perturbations, and 
finally, a stateme ,t of the overall objectives of this v/ork . 


1.1 Overview on the Method of Averages 

The purpose of the Method of Averages is to eliminate 
’’fast variables" from the differential equations. It is based 
on the following Main Theorem of Averaging (reference 1) for 
the periodic case: 

Theorem 

Consider the xnatial value problem 

-V -► 

— =ef(x,t,£) , x(0,e)=x (1) 

dt ° 

with t € [0,“>] , e e [O.E^] . X € G , an open bounded 

set in , and the following conditions are satisfied; 

(1) f is defined in a connected set G ; 

->■ 

(2) f is continuous and uniforraally bounded in G ; 





(3) f is Lipschitz-continuous with respect to x 

in g’*' ; 

(4) the li.nit 

lim ?(x,t,e) “ f(x,t,0) 
c-‘>-o 

exist,s uniformally in G ; 

(5) ? possesses a bounded derivative of e in G ; 
C6) ? is periodic in t with period T 

Consider the associated initial value problem 

dy ^ -»■ 

— = e ^o(y) , y(0,e) = x^ (2) 

dt ° 

where ^ 

f„(y) = - / f(y,i,0)dT 
T Jq 

If it is also assumed that y(t,e) exists, then 

5(t,e) - y(t,e) = QCe) 

1 

on the time scale — 

e 

Comments 

(1) Equation (1) is the "standard form" for the averaging 
method. Not all problems can be reduced to the standard 
form. 

(2) The conditions expressed above are generally satisfied 
by orbital mechanics problems. 


The first three conditions ensure that the solution to (1) 
exists and is unique. 
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(3) The theorem establishes that the solution of the Averaged 
Differential Equations (2) is always close to tho true 
solution, in the time interval of interest, 

(4) Notice that the right side of (2) contains only slowly 
varying variables. Therefore, these equations should be 
easier to solve. 

The Generalized Method of Averaging is discussed in ref- 
erence 2, and concerns the system of equations: 

- e f|Cx,y,e) , i = 1, 2, n 

dt ^ 


dy, ^ ^ ^ 

— ^ = w Cx) + e g. (x,y,F.) 
dt J 1 


j = 1, 2, • • ' , m 


where f, and g. are periodic in each component of y with 
^ J 

period 2 tt. The yj are called "rapidly rotating" phases. 

The purpose of the method is to eliminate y from the differ- 
ential equations. 

The Method of Averages has found widespread use in ap- 
plied mathematics and engineering. Different formulations, 
examples, as well as rigorous developments are given in 
references 3, 4 and 5. 

In this report, a Modified Method of Averages is devel- 
oped for application to equations of the form 

— = c f(x,y,t,e) 
dt 


dy 


= w(x) + e g(x, y, t ,e ) 


dt 
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Nota that tnere is only one rapidly rotating phase, but that 
the independent variable appears explicitly on the right hand 
side, It will be found that the KS equations can be put in 
the above form. For the applications presented in this report, 
it is convenient to eliminate only y from the differential 
equations. 


1,2 KS Total Energy Orbital Elements 

A very powerful method for the numerical solution of the 
differential equations of satellite motion is the total energy 
formulation of the KS element equations. The theoretical 
details of the KS method are developed in reference 6. Ref- 
erences 7 and 8 give evaluations of the KS elements when 
applied to numerical computation of satellite orbits. It has 
been found that the KS formulation offers the following ad- 
vantages : 

(1) Instabilities associated with solving the two-body 
(conic) equations are eliminated. 

(2) An orbital frequency based on the total energy 
gives more accuracy to calculations of the in- 
orbit (downrange) position, 

(3) The differential equations are "sm lothed" for ec- 
centric orbits because the eccentric anomaly is the 
independent variable. 

(4) The equations are less sensitive to roundoff and 
truncation errors in the numerical integration al- 
gorithm. 

It will be shown in Section 2 that the total energy fea- 
ture can be used to eliminate the tesseral geopotential terms 
from the numerical orbit computation. This is an additions! 
advantage of the KS formulation. 




where are the zonal terms and the tesseral terms. 

The zonal terms represent a potential function that is inde- 
pendent of the rotation of the earth. Thus, V.j. contains 

the longitude dependent part of the geopotential , and therefore, 
contains explicitly the time. 

It has been shown by analytical satellite theories that, 
whereas the zonal terms contribute long period and secular 
terms to the solution, the tesseral terms contribute primarily 
periodic terms. For solutions that are valid over several 
hundred revolutions, tesseral resonance terms must be included. 
However, for many applications, resonance effects are not im- 
portant and they will be neglected in this report. 

Even though the tesseral terms in the solution are small 
and periodic, it has been shown in references 9 and 10 that 
they can contribute large errors if neglected from numerical 
integrations. Downrange errors of 20 to 30 km are typical 
for near earth satellites (reference 9). These error.s can be 
attributed to the fact that an incorrect "mean" mean motion 
results when tesseral terms are arbitrarily dropped from the 
differential equations. 

It was shown in reference 10 that neglecting a small 
periodic perturbation from the differential equations, does 
not necessarily mean that the error in the solution will re- 
main small. Consider the following differential equations 
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dx 

■ — ■ “ Ea airtWt x( t»o) = x 

dt 


dy 

— “ X + eb ainuit y(t=*o) « y 

dt 


where e Is some small parameter. This situation is similar 
to the geopotential tesseral perturbations. The x is analo- 
gous to the mean motion and y is analogous to a fast vari- 
ble such as the mean anor.aly. Both are being perturbed by a 
periodic function. The solution for x and y is 


ta 

X =5 — (l-ocawt) + x„ 

0 

w 

c& f-.? eb 

y ** (x '*■ — ' t - — sintJt + — (l-ccsmt) + y 

0 .. ,2 .. o 


Suppose that the small perturbation is neglected. This 
is equivaxent to "crude" averaging, mentioned in reference 2. 
The solution is 


X ' “ X, 


y ' = X t + y 
^ 0 *^0 


The difference in the two solutions is 


x-x' * - (1-eosojt) 
a 


e ea eb 

v_y' _ t - — aintjit + — ( 1 - <?c>3£iJt ) 
2 

U1 UJ 


a 
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Note, that disregardlnK tho perturbation r(‘^;uIL.s tn a pcriodi<’ 
error In x . However, the error In y Is linear. 

Suppose, Instead, that an average value for x 


1 f ca. 

X’’ - — ; — (l-eoaO) •(• x^l dO 

2tt j ' tu " ’ 

o 


is inserted into the differential equation for y . Then, 
neglectinti the perturbation, the solution for y is 


or 


y" 




dt 


ca 

C — + 

ui 


* >’o 


The error incurred by this solution is periodic. 

This example and the numerical experiments of rt'feronce 9 
illustrate that crude averagiiiR will lead to incorrect results 
It will be shown in this report that with a prrjper initializa- 
tion of the mean motion, the numerical iniegration (uin jinx cc'd 
without includlnfr the tesseral geopotential terms. 


1.4 Qb.jectives of this VVork 

An important problem that is found when attempting the 
numerical solution of the satellite differential equations, 
is that the maximum size of the numerical Integration steps 
is limited by the high xrequency (short period) terms con- 
tained in the geopotential. This problem appears even with 
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the KS formulation. Even though their amplitudes are small 
and may he considered negligible for many applications, the 
short period terms cause the following practical problems: 

(i) Since the tesseral terms depend on time, the orbital 
frequency based on the total energy is no longer 
constant. If small steps are not talten to "frack" 
the high frequency oscillations, large down range 
errors can result. 

Hi) If the tesseral terms were simply neglected in the 
numerical integration, larger steps could be taken 
but unacceptable intrack errors would result 
(reference 9) . 

(Hi) Evaluation of all tesseral terms in the numerical 
integration force model consumes a ma,jor part of 
the computation time. 

An approach to solving the three above problems is dis- 
cussed in this report. Basically, the idea is to carry out 

t 

a numerical integration without the tesseral teims. In order 
to avoid the second of the above problems, the KS elements 
are initialized with a mean frequency (total energy) . The 
mean values are obtained via a specialized application of the 
Method of Averages. Since the force model contains only zonal 
geopotential terms plus, possibly, atmospheric drag and ex- 
ternal bodies, computer run time and stability problems asso- 
ciated with the tesseral terms can be avoided, 


t 


It will be assumed here that there are no important reso- 
nance effects. 
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2.0 MODIFIED METHOD OF AVERAGES FOR KS EQUATIONS 

One of the reasons for applying the Method of Averages 
to differential equations is to eliminate terms that do not 
make important contributions to the solution, but that can, 
nevertheless, cause difficulties in obtaining a numerical or 
analytical solution. That will be the objective in this sec- 
tion with regard to the KS differential equations. It was 
mentioned in Section 1.3 that neglecting the tesseral geopo- 
tential terms would cause a linear growth in the error of 
predicted satellite positions. In reference 10 it was shown 
experimentally that this linear error could be eliminated if 
a mean KS frequency could be found. In this section, a 
modified method of averages is developed and applied to the 
KS equations. This results in a mean KS frequency w that 
is defined by a partial differential equation. When this 
value of w is used in the numerical inte:*"ftion , the tesser- 
al geopotential terms can be safely neglectea for many appli- 
cations . 


2.1 

The 

KS 

Differential 

Equations 



The 

KS 

total energy 

element differential equations 

that 

will 

be 

used 

here are given in reference 6, section 23: 




dm 

r 3V 

1 




— 

ss _ — _ . 

— (u*, 

(3) 



dE 

3t 

8m 




di 

1 , 

^ 2 dm 



r -V . 4 - -i * 

I U - rV + 2(u, Q)J (u,u ) , (4) 

8o)^ dE 


dE 
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where 


da 


1 


2 dm 


E 

= 

— 

— 

Q " 

U 1 

sin 


dE 




m dE 

J 


2 

d^ 

- 


2 

dw 


E 



— 

— Q 

— 

— u 

oos — 


dE 

am® 

CO 

dE 


2 



L 









r 


1 





Q = - J l" 


3x 


(rV) 


+ - P 


( 5 ) 


( 6 ) 


(7) 


u , u , a , t and ^ are 4-vectors. 

$ is a 3-vector of perturbing accelerations that may or 
may not be derivable from a potential. V Is the disturbing 
potential function. 

For the purposes of this section, we will assume that 

->■ 

P is zero and that V is composed only of geopotential terms. 
In this case, equations (3) and (7) become 


dm r 3V 

— = ( 8 ) 

dE at 


^ = - i 


au 

The gravitational acceleration on a satellite is — 

, 

where 

U = - - V 
r 

and 

V = + V_ (10) 

^ 1 


— (rV) 
3x 


(9) 
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Thus the g©opc't®o^i0-l separated into zonal (time Inde- 

pendent) terms and l:essc ral (time dependent) terms, and 


R_ \ n 


P (ain p) 
n 


’^--y y P Uin ♦)fc oo. mx 

‘ V r/ L 


+ S sin ra X 
nm 


where 

4> - geocentric latitude, 

A = geographic longitude, 

= ordinars^ legendre function, 

P = associated legendre function, 
nm ’ 

= equatorial mean radius of the earth. 

C 

J , C , S = geopotential coefficients, 
n n m n m 

The time appears explicitly in V through A : 

T “ 

A = - (0^ t + (13) 

where 

ijj = inertial longitude of the satellite 
= Greenwich longitude. 

The time is given in terms of KS-elements by 

t=T--(u, u*) (14) 

Oi 
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Making use of equations (13), (i-l) and the transformation 
equations in section 19 of reference 0, it can be shown that 
may be v^-ltten as 



(15) 


( 16 ) 


and 



i-l 


Comiuents 

(1) Notice that has been brought outside the* summat ion . 

Therefore, the magnitudes of the funcrions A and B 

n m n in 

are of order 1, 

(2) The functions A and B are periodic in E with 

nir nm 

period 2tt 

(3) V^, is "purely" periodic in 0 with period 2tt . That 

means that there are no terms in which are indepen- 

dent of 0 

Making use of the expressions for V . V , assuming 

z T 

that ? = 0 , and using the new time element O . the KS 
differential equations can be formally expressed as 
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— = e" (u),0.a.|,E) (17) 

dE 


d0 u . ^ ^ 

— » 6 + e g, (u,a,S,E) + e® (u , 0 , a , 6 , E) , (18) 

dE Sw® *' 


— = e h, (w,a,3,E) + £* h„ (aj,G,a,B,E) , (19) 

dE ^ 


e (u,a,f,E) + (w,0,a,^,E) . (20) 

“V 

The parameter e = J 2 has been introduced. The function 
and kj^ depend only on and the functions h 2 and k 2 

depend only on V.^, . Based on the above comments concerning 

, it can be concluded that the functions f 2 . ’ 

112 and 1^2 are purely periodic in Q , i.e. 



2 TT 


f2 (w,0,a,|,E) d 0 = 0 


( 21 ) 


2it 


g2 (ai,0,a,^,E) d 0 = 0 




^2 (w,0,a,l,E) d 0 = S 


( 22 ) 


(23) 


0 
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1*2 (w.e.aJ.E) d 0 « 5 (24) 

These i-esults will be important in the theoretical develop- 
ments with the Method of Averages. 

2.2 The Modified Method of Averages 

For the developments that follow, there will be no loss 
of generality if the order of the KS differential equations 
(17) through (2'"’) is reduced. Thus, the equations will be 
represented as 

dw 

— * f- (w,0,a,E) (25a) 

dE 


dG 

— = Li(w) + e g, (oj.a.E) + g„ (u,0,a,E) (25h) 

dE ^ ^ 


da 

— = e h (tjj,a,E) + c* h„ (w,0,a,E) (25c) 

dE ^ ^ 

where 

u(ai) = 0_ (26) 

So)3 E 

The functions on the right side of equations (25) have the 
properties; 

(i) gj , g 2 , hj, li 2 are periodic functions of E 

Cii) g 2 and li 2 are purely periodic functions of 

0 . They represent the tessera 1 terms of the 

geopotent ia 1 . 
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(iii) gj^ and do not depend on Q . They represent 

the zonal terms of the geopotential. 

Since the right sides of equatlt.ns (26) are periodic in 
the angular variables, it is appropriate to apply averaging 
for their solution. A Modified Method of Averages (MMOA) 
will be applied in this section to equations (25). The 
MMOA is a modification of the averaging method that is de- 
veloped in reference 2. 

The purpose of the MMOA is to eliminate only those 
periodic terms that depend on a depandsnt angular variable 
(i.e. time). This approach will eliminate the tesseral geo- 
potential terms from the problem. The averaged differential 
equations will still depend on the Independent variable E 
and the right sides will contain short period terms. However, 
numerical solution of the averaged equations will be enhanced 
because they will have the following properties: 

(i) The frequency of the fast variable will be a con- 
stant. This means that the numerical solution will 
have increased stability, i.e. larger steps can be 
taken . 

(ii) Evaluation of the right sides will be much quicker 
since the time consuming tesseral terms need not 
be computed. 

Given the differential equations (25), assume asymptotic 
expansions of the form 

(It = ( 1 ) + (w,0,a,E) + e^Ti 2 (w.Q.ct,^) + **' , (27a) 

0 = 0+ (IjJ,0,a,E) + (w,0,a,E) + ••• , (27b) 

a =a + exj (w,0,a.E) + c^X 2 (w,0.a.E) + ••• . (27c) 

The functions x^ . (i = 1,2, are required to be 

periodic in 0 and E with period 2tt . The new (averaged) 
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elements are defined by the differential equations, 


dw 

— “ eP, (m,a,E) + e^P, (w,a,E) + ••• , (28a) 

dE ^ 


dS" _ __ _ _ 

= UCw) + eG, (w,a,E) + e^G„ (w,a,E) + ••• , (28b) 

dE ^ ^ 


da 

— = eH, (w,a,E) + e^H„ (w,a,E) + ••• (28c) 

dE ^ ^ 

Equations (28) are referred to as the "averaged" equations. 

It will be shown that they have the properties that were pre- 
viously discussed. 

The functions (i « 1,2, •••) 

are not yet known and will now be determined by making expan- 
sions in powers of e . The coefficients of the powers of 
e will provide a set of partial differential equations. The 
solution of these equations will give the required functions. 

Power series expansions in e are obtained by the fol- 
lowing algorithm: 

(i) Compute the total derivative of both sides of equa- 
tions (27), using the chain rule for the functions 
on the right side. 

(ii) Insert equations (28) into the right sides of the 
equations produced in step H) . This will give 


dw 

d0 

da 

expressions for — 


, and — in terms the 

dE 

dE 

dE 

averaged variables. 



Insert equations (27) 

into 

equation (25) and expand 


in powers of e 
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(iv) Equate the expressions for tho total derivatives of 
0 ), 0, and a that were produced in step fti) to 
those that were produced in step (Hi). 

(v) Collect terms of like powers of c . The coeffi- 
cients are the required partial differential equa- 
tions . 

The defining equations are: 

£ ° Terms 

U(w) = u(u) (20) 

e ^ Terms 

3n, 3n, 

+ U — - + — ~ = 0 C30) 

^ 90- 

9$, 9t(i, 9u 

G + U — - + — ~ = n, — + g, (31) 

^ 95^ 9E ^ 3w ^ 



9n, 9n, Sho 

F + F + G + H, + U ^3“ + — - 

^ ^ 9(1) He 9a ae 9 e 


(32) 


( 33 ) 
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94>, 3<|)i S'f'i 

G, + F, ^ + G, + H, “3^ + U •'■ — ~ 

^ ^ 3w 30 3a 30 3E 


(34) 




3u 1 3*U 1 3k 

n, — +-nf — +_n,— + “ + «2 


3w 2 ^ 3u* 2 '■3w 


3a 


E>Xi 

H + F — i- + G. 
^3(0 


9Xi 

30 


+ U 




1 3h 1 3h 

" “ 1 1 ~ir ■*' - X 1 -zr 

2 3u) 2 3a 


(35) 


Note that all dependent variables in equations (29) through 
(35) are barred. 

Equation (30) has two unknown functions and n^ 

Thus, define 

= 0 (36) 

Then, ri| can also be set to zero, i.e, 

= 0 (37) 

In equation (31), make use of equation (37) and define 


G 


1 


1 

2tt 


.2it 




d 0 


'0 
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Since doss not depend on G 

Gj (w.a.E) « g^Cu.a.E) (38) 

The defining equation for <|>j is therefore 

3 (|) 34 > 

U — » 0 
30 3E 

30 'that 

“ 0 (30) 

Similarly, from equation (32), 

II^ (w,ct,E) = hj^ (w,a,E) (40) 

and 

Xj = 0 l41) 

There has been established the following important results: 
(i) All first order terms in the transformation equa- 
tions (27) are zero. 

(ii) The averaged KS frequency w is constant to 
first order. 

(Hi) The first order terms in the derivatives of the 
remaining averaged KS elements are unchanged 
from the original equations. 

Now the second order equations (33), (34) and (35) will 
be used. Making use of “ 0 , equation (33) becomes 
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Sn, 

— £ + — £ M 


F + U 

^ 3^ 3E 


The above equation has two unknown fur tlons, ao one can be 
arbitrarily defined. Remember requir( i to be purely 

periodic in S’ and/or E . Therefore, set 


2iT„2ir 


F, 





dSdE 


0 •'0 


But recall that is purely periodic in S (since it con- 
tains only tesseral terms). Thus, is zero and 02 is 

defined by the equation 


3n, 3n, 

U + — £ = 

30 3E 


(42) 


Making use of similar arguments, it is found that 


02 - 0 


(43) 


II 2 = 0 


(44) 


and 


^2 J^nd X 2 defined, respectively, by the equations 


3^2 ^‘^*2 • 

30 3E ^ 3w 


+ g. 


(45) 


U 


3 X 2 . ^^2 
+ = h. 


(46) 


30 3E 
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Th® additional important results have been established: 

(i) The derivative of w is zero through second order 
in E 

(ii) Tlie averaged differential equations do not depend 
on the tesseral terms, through second order In t 
(Hi) Transformation between osculating and mean elements 
is givei by the solution of partial differential 
equations . 

If the power series expansion method is carried to third 
order, the following results involving w are obtained: 

F 3 = 0 , (47) 


(48) 


Thus, [D is constant through third order in t , The func- 
tion is defined by a partial differential equation, the 

right hand side of which contains coupling effects between 
tesseral and zonal geopotential terms. 

Finally, the averaged differential equations are: 
dw 

— = OCeM (49a) 

dE 


dG _ 

— = U(oj) teg, (w,a,E) + 0 (e®) , (49b) 

dE ^ 


da __ _ 

— = e h (w,a,E) + 0 (e:^ ) 
dE ^ 


3n, 


3p, 


U 


3 _ 


3G 


3E 


G, 


9n, 


ao 


+ H 


3ri2 
^ 3a 


(49c) 
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The transformation equations betwoon mean and instantaneous 
elements are: 

03 “ w + E* U 2 <w,S',a,E) + (50a) 

0 = 0 + E^ (|)2 (w,8,a,E) + 0(e®) » (50b) 

a = cT + Xj (m.Oja.E) + Q(^^) (60c) 

Based on the averaging theorem discussed in Section 1.1, 

it is expected that the averp.ged solution will differ from 
the exact solution by order e 


)utational 


It will be necessary to express the averaged elements in 
terms of the instantaneous elements. The required transform- 
ation can be obtained from equations (50) by observing that 


(w,0,a,E) = Ti2 (f^,^,ot,E) + O(e^) 


with similar expressions for and X 2 - Therefore, the 

required transformation equations are obtained by "reversing" 
(60): 

w = u3 + E^ ri 2 (t^ ’.a,E) + Q(e^) i (51a) 


0 = 0 + E^ ())2 (o),0,a,E) + , (51b) 


a = a + X 2 (w,0,a,E) + 0 (e^ 


(51c) 
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The computational algoritam :Por solving the differential 
equations is as follows; 

(X) Solve the partial differential equations (42), (45) 
and (46), Notice that the equations depend only on 
the tesseral geopotential terms and that the oscu- 
lating values of the elements w , 0 , a are 

used. This is allowed because of equations (51), 

(2) Use the values of and X 2 obtained in 

step (1) to obtain w , 0 and a by equations 

(51). 

(3) Solve the averaged differential equations (49) by 
using the initial conditions obtained from step (2), 
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3.0 APPLICATIONS OF AVERAGED KS ELEMENTS 

In this section, the averaged KS equations will be 
applied to satellite orbit prediction problems, By making 
use of the theory developed in Section 2, it will be found 
that the objectives stated in Section 1.4 are realized. 


3.1 Numerical Solution of the Partial Differential Equation 

It has been shown that the averaged elements were deter- 
mined by functions that are solutions of partial differential 
equations. Nothing has been said, so far, about the solution 
of these equations. A numerical solution algorithm will be 
developed in this section. 

First, an important simplification can be made. Consider 
the averaged differential equations (49) and make the follow- 
ing observations: 

(i) The right sides of equations (49) do not depend on 
0 . Therefore, any small initial error in 0 

vi.e. in time) will remain constant and will not 
affect the solution, 

(%i) If a is replaced by a in and h^^ , the 

error is of order . Since this error is of the 

same order as neglected terms, in equations (49b) and 
(49c), it is permissible to insert a into and 

h^ . 

(iii) The function U is a secular term depending only 

on w . Therefore it is important that w be ac- 
curate, otherwise there would be a secular (linearly 
increasing) error. This is the source of the error 
demonstrated in reference 9. 

(iy) Since w is a constant, it is necessary only to 
initialize lo , using equation (51a). 
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Protn the above considerations, it is concluded that only the 
partial differential equation (42) for ri 2 needs to be solved. 

Consider the equation for r\^ : 

3rin ^ . 

U ^ = f (u,0,a.1,E) (52) 

3G 3E 


where the bars have been removed from the dependent variables 
since, according to Section 2.3 and equation (51a), it is al- 
lowed to express in terms of the unaveraged elements. 

Also, the dependence of f^ on all the KS elements is shown 

in equation (52). 

A method of solving equation (52) is to express the right 
side as a truncated double fourier series in the variables G 
and E , This idea is presented in section 28 of reference 6 
for a similar problem. The expression for f 2 is then 


1 1 ^ 

P (0,E)2 - a + - 2-/ (a,- 

2 4 0° 2 ^ 


eos i 0 + a. sin i 0) 
o lO 


M 

- . 

O OJ 


2 J-1 


00 a j E + b . a in j E) 

oj 


L M 


; EE 


(a,, - b..) COB (iO+JE) 
ij ij 


i=l j = l *- 

+ (a. . + b ) aoa (x0-jE) 

ij ij 

+ (b.. + a..) ain (10+jE) 

ij ij 


( 53 ) 


(b,. - a..) sin (iO-jE) 
13 X3 
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The coefficients a. 


a, 


^ii ' "ij ' ’ "ij 

^ and w only. Numerical algorithms for computing the co- 
efficients are given in the Appendix. These algorithms make 
use of the fact that f_ is a known function and can be com- 


depend on a 


puted by the routines in the KSFAST program (reference 11). 
See also the discussion given on page 163 of reference 6. 


It is known that does not have a constant term or 

any terms independent of 0 . Therefore, some of the coef- 

ficients must be zero, i.e. 


a ** 0 (54) 

oo 



j = 1,2,3. 


The solution of equation (52) is then 

L 

1 1 _ 
n- (0,E) = — / j — (a, sin i 0 - a, cos i 0) 

^ 2U 7:7 i 




i=i j = i 


a. -b, . 
(iU+j) 


stn 


(iO+jE) 


( a . . +b . . ) 

+ — U .- , IX - sin (iQ-jE) (55) 

(iU-j) 


y ^ .f » / 

^ 00 s (10+jE) 

iU+.i 


(b. .-a. . ) 

+ — cos (i0-jE) 
iU-a 
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3.2 Computat lonai A ikorithm 

According to the discussion in Section 3.1, it i orly 
necessary to compute the initial value of w ; 

where is the Initial KS frequency computed from the 

initial position, velocity and time. It is assumed that 

E = 0 . Also, 

o 


For a complete discussion of the initialization of the KS 
elements, refer to section 19 of reference 6. 

Inserting 0=0^ and E = 0 into equation (55), the 

required expression is 


1 1 

2 0 2U fr* 1 ^ 

1-1 


sin 1 Q -a, ooai0) 
0 io o' 


I M 


(5b) 


(jb^j+iUa^j) ain i 0^ 


+ (Jb^j-iUa^j) 008 i 0Q 


the 


The compucation of satellite orbits is carried out 
following steps: 

(1) Compute the initial KS elements 


in 


*^0 
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( 2 ) 

(3) 


from the initial satellite state vector, by mak- 
ing use of the transformation equations (reference 6, 
ection 19). Be sure that the disturbing potential, 
function that is used to compute includes 

all the tesseral and zonal terms in the geopi Lential. 
Compute the coefficients a^^ , b^^ 

using the equations in the Appendix. 

Evaluate equation (58) for 


■ '‘u ■ ^13 


and then com- 


pute 


- 03 ^ - n2<0o.O) 


(4) Replace m with w in the set of initial KS 
elements . 

(5) Carry out the numerical integration of the KS dif- 
ferential equations (3), (4), (5), (6) with the con- 
dition that the function V contains only the zonal 
geopotential terms. Atmospheric drag and third body 
perturbing accelerations are included in 5 

A computed trajectory based on the above algorithm will 
agree to within a few hundred meters to the true trajectory. 


3.3 Satellite Trajectory Prediction Experiments 

A set of numerical experiments has been carried out using 
the averaging option (tesseral initialization) in the KSFAST 
program (reference 11, section 4.9). Selection of this option 
causes execution of the satellite orbit prediction algorithm 
that was described in Section 3.2, The experiments described 
in this section demonstrate the computation time savings that 
are available with this option. 

Two sets of initial conditions are considered. Perturb- 
ing effects on the satellite orbit are atmospheric drag and 
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a nonspher ical , nonhomogenous earth. The genpotent ial nottii- 
tmlly includes all terms through eighth ordc?r and eighth 
degree. The purpose of the comparisions is to demonstrate 
that if the averaging algorithm is used, computation time can 
be reduced by a factor of about 4 . This is because the 

tesseral geopotential terms no longer need to bo included in 
the numerical Integration, 

The initial conditions for the two trajectories are: 


Case 1 

Altitude 
Eccentricity 
Incl inat ion 
Longitude 
Latitude 

Epoch January 1, 1977 

Case 2 

Perigee Altitude 

Apogee Altitude 

Eccentricity 

Inclination 

Longitude 

Latitude 

Epoch January 1, 1977 


300 km 
0 

28 deg. 
136.93 deg. 
0 

22 hrs. 


300 km 
500 km 
.015 
28 deg. 
136.93 deg, 
0 

22 hrs. 


Figures 1 and 2 .show how the error in position will in- 
crease if the tesseral terms are not included in the numer- 
ical integration and no averaging is done. (Compare with the 
experiments in reference 9.) This error is almost entirely 
in the along-track direction and is caused by an incorrect 
"mean" mean moti n. Therefore, the error grows linearly, 

The purpose of the averaging method is to remove this linear 
error, A periodic error will still remain, but will be small 
(=200 m. ) . 






FIGURE 1: Position Error (Ko Avora^o), Cnsa 1 



0 4 8 12 IQ 20 24 SS 32 

Humber of Revolutions 

FIGURE 2: Position Error (Ho Averngo), Csso 2 
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Results of satellite trajectory prediction experiments 
are displayed in Tables I and II. Two types of comparisonB 
are made; 

(1) Position error - The errors for the averaged method 
are shown, as well as the No Average case. In both 
cases, tesseral terras are neglected. However, in 
the No Average case, the averaging is not done. 

Almost all of the position error in each case is In 
the along-track direction. 

(2) Computation time - The times are shown for the 
Averaged case (no tesseral terms) and a Precl.sion 
case (all tesseral terms included in the numerical 
integral ion ) . 

Discussion of Results 

(1) The averaging method effectively reduces the linear 
growth of the along-track error in position. Periodic 
errors remain and are on the order of a few hundred 
meters . 

(2) The larger errors found in Case 1 occur because of 
coupling between drag and tesseral terms. Case 1 
is a rather ' ow circular orbit and the coupling ef- 
fects are important. For Case 2, the satellite 
spends less time in the dense atmosphere so the 
coupling effect is not so important. Even so, the 
coupling is only important after about four days. 

(3) The errors in the Averaged Method in both cases are 
so small that they would not show up on the scale 
of Figures 1 and 2. 

(4) The run time comparisons of Tables I and II show 
that the Averaged Method is about three to four 
times more efficient than the Precision Method. 

(5) About two seconds are required for the averaging 
initialization algorithm. Therefore, the Averaged 
Method is most efficient for prediction intervals 
of more than one day. 
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TABLE I.- NUMERICAL COMPARISONS, CASE 1 


Prediction 

Interval 

(Days) 

Position Error (km) 

Computation Time (sec) 

Averaged 

Method 

No 

Average 

Averaged 

Method 

Precision 

Method 

1 

.10 

21.4 

3.5 

8.8 

2 

.27 

42.9 

5.0 

17.7 

4 

,80 

86.3 

9.9 

35.4 

0 

1 .58 

130.3 

14.1 

53.1 

8 

2.74 

174.8 

18,4 

70.7 

10 

- 4 .40 

220.0 

22,6 

88 .4 

15 

11.00 

337.0 

33.2 

132.6 


TABLE II.- NUMERICAL COMPARISONS. CASE 2 


Prediction 
In terval 

Position 

Error (km) 

Computation 

Time (sec) 

(Days) 

Averaged 

Method 

No 

Average 

Averaged 

Method 

Prec i si on 
iiethod 

1 

.15 

21.2 

4.0 

8.0 

2 

.29 

43.2 

6.1 

17.3 

4 

.63 

85.5 

10.3 

34.6 

1 

0 

.41 

126.7 

14 .4 

51.9 

8 

.29 

169.1 

18.6 

69.1 

10 

.48 

214 .4 

22.7 

CD 

00 

15 

.93 

' 319.2 

33.1 

129.6 
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(6) The efficiency of the Averaged Method is because 
the tesseral terms do not need to be included in 
the step by step numerical integration. Also, a 
large stepsize can be taken because the averaged 
differential equations are more stable. This is 
because the KS frequency element io Is nearly 
constant . 

The errors of a few hundred meters in the solution by 
the Averaged Method are caused by the diurnal (dally) varia- 
tions produced by the neglected tesseral terms. It would be 
possible to add these variations to the solution by convert- 
ing back to osculating elements via equations (51), The 
position error would then be reduced to a few meters. This 
procedure could be added to the Averaged Method in case more 
accuracy is required. There would be a negligible increase 
in computation time since the back conversion would be done 
only for output. 

A further refinement would eliminate much of the long 
term error growth of the Case 1 Averaged Method. As was men- 
tioned above, this error is due to coupling between atmospheric 
drag and tesseral perturbations In the altitude. Since the 
density of the atmosphere increases exponentially with decreas- 
ing altitude, the small variations in altitude caused by tes- 
seral perturbations become important over many revolutions. 
Equations (51) could be used to compute the actual altitude 
at each numerical integration step. This altitude would go 
into the density model. More accurate calculation of drag 
perturbations would I'esult . 
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4.0 CONCLD'SIDNS 


The Averaging Inlt ializinLlon Algorithm that is rievol op,«d 
In this report can clecreaHR the conputation time required for 
predicting near earth satellite orbits. For prediction intt^r- 
vals of one or more days, the computation time will be reduced 
by a factor of three to four. The? methov., is particularly 
suited to iterative applications where accurate in-orblt sat- 
ellite position Information is required. 
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APPENDIX 

TRIGONOMETRIC POLYNOMIAL APPROXIMATION IN TWO VARIABLES 

Suppose that there exists a function f(0,(()) which is 
2n periodic in both variables. It is required to establish 
a fourier approximation for f in the region 

O<0<2ir '0<4><2ti 

The approximation will be based on numerical evaluations of 
f at an odd number of values of the arguments. Therefore, 
let 


x=0, 1, 2, •••,2L 


y = 0 , 1 , 2 , • • * , 2M 


where 


2L+1 

X = 0 

2'tt 


2M+1 
y = 

2 TT 


(Al) 


(A2) 


The fourier approximation formula in one dimension is 
(reference 12) 

L 

1 2v \ 

f(x) ~ — A + / y | A. cos kx B, sin kx) , (A3) 

2 ” iTTlV I* 2L+1 2L+1 / 


where the coefficients are 


A. 

.1 


21. 


2L+1 x=o 


f(x) 


2 71 

cos j 

2L+1 


X 


j = 0,1,..., L, (A4) 


PRECEDia<3iiPAGaLBLAHIC*taOT*|iLM£D. 
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2L 


B. = 


J 2L-M ^ 


X)f(x) 


etn 


2ir 


2L+1 


jx , j»l,2,*‘*,L. (A5) 


Assume that ihe coefficients A, and. B, in (A3) are functions 

k k 

of y . Then they car also be expressed as fourier approxi- 
mations: 


1 

Aj(y) s - Aj., + 


K / 

s(' 


2 IT 
2M+1 


£.'0.9 j y 


2it 

+ b, , sin J 


Ij 


2M+1 


i . 


i = 0, 1, 2,**‘‘L, 


(A6) 


M 

B,(y) - - a.^ + . cos 


1 

2 


.1 i 


2n 


i'M+1 


J y 


(A7) 


2tt 


+ b. . sin J 

2M+1 


i y^ . i = 1. 2, 


, L. 


Expressions (A6) and (A7) are inserted into (A3). Using equa- 
tion (A2) and trigonometric identities, the double fourier 
approximation is 

L 


1 1 

i(0»<f>)2~a +~/j(3-, cosiO + a. nin i Q) 

oo o io lo 

4 2 i=] 


+ 


M 

- (a . 

mmmid O J 


2 Pi 


cosj0 + b , stnj(()) 
o.) 


(AS) 


L M 


EL 

2 i=l j=l 


- b.^ ) cos (i 0 + j (j)) 





